function compareCtdsMainBasin(modeldir,modelbasename,year,outfile,...
							  do_loadObs,do_loadModel,do_graphics);
							  
% compareCtdsMainBasin(modeldir,modelbasename,year,outfile,...
%					   do_loadObs,do_loadModel,do_graphics);
%
% neil banas, 2014

Tdir = toolstart;
ctdDir = [Tdir.data '/obs_data/ctd_bottle/'];

if nargin<1, modeldir = '/pmraid3/sarahgid/runs/salish_2006_4/OUT_lp/'; end
if nargin<2, modelbasename = 'ocean_his_lp_'; end
if nargin<3, year = 2006; end
if nargin<4, outfile = [Tdir.output 'compareCtdsMainBasin']; end
if nargin<5, do_loadObs = 1; end
if nargin<6, do_loadModel = 1; end
if nargin<7, do_graphics = 1; end

Tcontours = 0:20;
TcolormapLevels = [0 6 18 20];
Scontours = [0:4:28 28.25:0.25:34];
ScolormapLevels = [0 28 32 34];
dTcontours = 0:0.5:6;
dScontours = 0:0.5:16;
dTdScolormapLevels = [0 0.05 0.5 1, 0.01];


if do_loadObs
	% load observations (T,S) along the main channel
	load stdTracksAndStns jdf2psTrack;
	[ex,fnames] = obs_extract(ctdDir,...
		{'salinity','temperature'}, ...
		[datenum(year,1,1) datenum(year+1,1,1)],...
		'section', jdf2psTrack, 5);
	casts = obs_separateCasts(ex);
	for i=1:length(casts)
		mb.x(i) = casts(i).x(1);
		mb.y(i) = casts(i).y(1);
		mb.t(i) = casts(i).t(1);
		mb.dist(i) = casts(i).dist(1);
		f = find(~isnan(casts(i).temperature+casts(i).salinity));
		fsurf = f(abs(casts(i).z(f)) == min(abs(casts(i).z(f)))); 
		fbot = f(abs(casts(i).z(f)) == max(abs(casts(i).z(f)))); 
		if ~isempty(f)
			mb.T0(i) = casts(i).temperature(fsurf);
			mb.Tbot(i) = casts(i).temperature(fbot);
			mb.S0(i) = casts(i).salinity(fsurf);
			mb.Sbot(i) = casts(i).salinity(fbot);
		else
			mb.T0(i) = nan;
			mb.Tbot(i) = nan;
			mb.S0(i) = nan;
			mb.Sbot(i) = nan;
		end
	end
	save(outfile,'mb');
else
	load(outfile,'mb');
end


if do_loadObs & ~do_loadModel
	warning('overwrote output file, have to load model again');
	do_loadModel = 1;
end
if do_loadModel
	% load model fields along the main channel
	ser = roms_createSeriesDef(modeldir,modelbasename);
	[T,mmb] = roms_extract(ser,'temp',ser.nctime,'profile',jdf2psTrack.y,jdf2psTrack.x);
	mmb.dist = repmat(...
				  reshape(jdf2psTrack.along_dist,[1 1 length(jdf2psTrack.along_dist)]),...
				  [size(mth.xm,1) size(mth.xm,2) 1]);
	mmb.T = T;
	mmb.S = roms_extract(ser,'salt',ser.nctime,'profile',jdf2psTrack.y,jdf2psTrack.x);
	save('-append',outfile,'mmb');
else
	load(outfile,'mmb');
end



if do_graphics
	[year,~] = datevec(mb.t);
	year = year(1);
	
	tm = squeeze(mmb.tm(:,1,:));
	distm = squeeze(mmb.dist(:,1,:));
	
	figure
	orient landscape
	
	subplot 221
	scatter(mb.t,mb.dist,40,mb.T0,'filled');
	caxis(Tcontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('surface T');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	subplot 222
	scatter(mb.t,mb.dist,40,mb.Tbot,'filled');
	caxis(Tcontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('bottom T');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	subplot 223
	contourf(tm,distm,squeeze(mmb.T(:,end-1,:)),...
		Tcontours,'edgecolor','none');
	caxis(Tcontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('surface T');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;

	subplot 224
	contourf(tm,distm,squeeze(mmb.T(:,2,:)),...
		Tcontours,'edgecolor','none');
	caxis(Tcontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('bottom T');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	colormap(tailedrainbow(TcolormapLevels));
	print('-dpdf',[outfile '_T']);
	
	figure
	orient landscape
	
	subplot 221
	scatter(mb.t,mb.dist,40,mb.S0,'filled');
	caxis(Scontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('surface S');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	subplot 222
	scatter(mb.t,mb.dist,40,mb.Sbot,'filled');
	caxis(Scontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('bottom S');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	subplot 223
	contourf(tm,distm,squeeze(mmb.S(:,end-1,:)),...
		Scontours,'edgecolor','none');
	caxis(Scontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('surface S');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;

	subplot 224
	contourf(tm,distm,squeeze(mmb.S(:,2,:)),...
		Scontours,'edgecolor','none');
	caxis(Scontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('bottom S');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	colormap(tailedrainbow(ScolormapLevels));
	print('-dpdf',[outfile '_S']);

	figure
	orient landscape
	
	subplot 221
	scatter(mb.t,mb.dist,40,mb.T0-mb.Tbot,'filled');
	caxis(dTcontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('T stratification');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	subplot 222
	scatter(mb.t,mb.dist,40,mb.Sbot-mb.S0,'filled');
	caxis(dScontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('S stratification');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	subplot 223
	contourf(tm,distm,squeeze(mmb.T(:,end-1,:)-mmb.T(:,2,:)),...
		dTcontours,'edgecolor','none');
	caxis(dTcontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('T stratification');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;

	subplot 224
	contourf(tm,distm,squeeze(mmb.S(:,2,:)-mmb.S(:,end-1,:)),...
		dScontours,'edgecolor','none');
	caxis(dScontours([1 end]));
	axis([datenum(year,1,1) datenum(year+1,1,1) -20 130]);
	set(gca,'ydir','reverse');
	title('S stratification');
	ylabel('Distance from AI (km)');
	datetick('x','mmm yy','keeplimits');
	colorbar;
	
	colormap(tailedrainbow(dTdScolormapLevels));
	print('-dpdf',[outfile '_dTdS']);
end